#include "stdafx.h"
#include "hnum_slucmn.h"

namespace harlinn
{
    namespace numerics
    {
        namespace SuperLU
        {

            void
            pxgstrf_super_bnd_dfs(
		                  const int  pnum, /* process number */
		                  const int  m,    /* number of rows in the matrix */
		                  const int  n,    /* number of columns in the matrix */
		                  const int  jcol, /* first column of the H-supernode */
		                  const int  w,    /* size of the H-supernode */
		                  SuperMatrix *A,  /* original matrix */
		                  int        *perm_r,   /* in */
		                  int        *iperm_r,  /* in; inverse of perm_r */
		                  int        *xprune,   /* in */
		                  int        *ispruned, /* in */
		                  int        *marker,   /* modified */
		                  int        *parent,   /* working array */
		                  int        *xplore,   /* working array */
		                  GlobalLUBase_t *Glu /* modified */	      
		                  )
            {
            /*
             * -- SuperLU MT routine (version 1.0) --
             * Univ. of California Berkeley, Xerox Palo Alto Research Center,
             * and Lawrence Berkeley National Lab.
             * August 15, 1997
             *
             * Purpose
             * =======
             *
             * Performs a symbolic structure prediction on a supernode in the Householder
             * matrix H, with jcol being the leading column.
             *
             */
                //GlobalLUBase_t *Glu = pxgstrf_shared->Glu; /* modified */	      
                register int krep, chperm, chrep, kchild;
                register int invp_rep; /* "krep" numbered in the original A */
                register int krow, kperm, xdfs, maxdfs, kpar;
                register int fsupc, k, jj, found;
                register int nrow;  /* union of the nonzero rows in a supernode */
                NCPformat    *Astore;
                int          *asub, *xa_begin, *xa_end;
                int          *xsup, *xsup_end, *supno, *lsub, *xlsub, *xlsub_end;

                /* Initialize pointers */
                xsup       = Glu->xsup;
                xsup_end   = Glu->xsup_end;
                supno      = Glu->supno;
                lsub       = Glu->lsub;
                xlsub      = Glu->xlsub;
                xlsub_end  = Glu->xlsub_end;
                Astore   = (NCPformat*)A->Store;
                asub     = Astore->rowind;
                xa_begin = Astore->colbeg;
                xa_end   = Astore->colend;
                nrow = 0;
                found = n + jcol;

                /* For each column in the H-supernode */
                for (jj = jcol; jj < jcol + w; ++jj) {

                  /* For each nonz in A[*,jj] do dfs */
                  for (k = xa_begin[jj]; k < xa_end[jj]; ++k) {
	              krow = asub[k];
	  
	              /* krow was visited before, go to the next nonzero. */
	              if ( marker[krow] == found ) continue;
	  
  	              /* For each unmarked nbr krow of jj ...   */
	              kperm = perm_r[krow];
	  
	              if ( kperm == EMPTY ) { /* krow is in L */
	                  marker[krow] = found;
	                  ++nrow;
	              } else {
	                  /* krow is in U: if its supernode-rep krep has been explored,
		                           skip the search.       */
	                  krep = SUPER_REP( supno[kperm] );
	                  invp_rep = iperm_r[krep];
		  
	                  /* Perform dfs starting at krep */
                          if ( marker[invp_rep] != found ) {
		              marker[invp_rep] = found;
		              parent[krep] = EMPTY;
		              if ( ispruned[krep] ) {
		                  if ( SINGLETON( supno[krep] ) )
			              xdfs = xlsub_end[krep];
		                  else xdfs = xlsub[krep];
		                  maxdfs = xprune[krep];
		              } else {
		                  fsupc = SUPER_FSUPC( supno[krep] );
		                  xdfs = xlsub[fsupc] + krep-fsupc+1;
		                  maxdfs = xlsub_end[fsupc];
		              }

		              do {
		                  /* For each unmarked kchild of krep ... */
		                  while ( xdfs < maxdfs ) {
			              kchild = lsub[xdfs];
			              xdfs++;
			              if (marker[kchild] != found) { /* Not reached yet */
			                  chperm = perm_r[kchild];
			      
			                  if ( chperm == EMPTY ) { /* kchild is in L */
				              marker[kchild] = found;
				              ++nrow;
			                  } else {
				              /* kchild is in U: 
				               *   chrep = its supernode-rep. If its rep
				               *   has been explored, skip the search.
				               */
				              chrep = SUPER_REP( supno[chperm] );
				              invp_rep = iperm_r[chrep];
				  
				              /* Continue dfs at snode-rep of kchild */
				              if ( marker[invp_rep] != found ) {
				                  marker[invp_rep] = found;
				                  xplore[krep] = xdfs;
				                  xplore[m + krep] = maxdfs;
				                  parent[chrep] = krep;
				                  krep = chrep;/* Go deeper down G(L^t) */
				                  xdfs = xlsub[krep];     
				                  maxdfs = xprune[krep];
				                  if ( ispruned[krep] ) {
					              if ( SINGLETON( supno[krep] ) )
					                  xdfs = xlsub_end[krep];
					              else xdfs = xlsub[krep];
					              maxdfs = xprune[krep];
				                  } else {
					              fsupc = SUPER_FSUPC(supno[krep]);
					              xdfs = xlsub[fsupc] + krep-fsupc+1;
					              maxdfs = xlsub_end[fsupc];
				                  }
				              } /* if */
			                  } /* else */
			              } /* if... */
		                  } /* while xdfs < maxdfs */
		      
		                  /* krow has no more unexplored nbrs:
		                   *    Place snode-rep krep in postorder dfs, if this 
		                   *    segment is seen for the first time. Note that
		                   *    the "repfnz[krep]" may change later.
		                   *    Backtrack dfs to its parent.
		                   */
		                  kpar = parent[krep]; /* Pop stack, mimic recursion */
		                  if ( kpar == EMPTY ) break; /* dfs done */
		                  krep = kpar;
		                  xdfs = xplore[krep];
		                  maxdfs = xplore[m+krep];
		              } while ( kpar != EMPTY ); /* do-while - until empty stack */
	                  } /* if */
	              } /* else */
                  } /* for each nonz in A[*,jj] */
                } /* for jj ... */

                DynamicSetMap(pnum, jcol, nrow * w, Glu);
    
            /*    for (i = 1; i < w; ++i) Glu->map_in_sup[jcol + i] = -i;*/
    
            #if ( DEBUGlevel>=1 )
                printf("(%d) pxgstrf_super_bnd_dfs(): jcol= %d, w= %d, nrow= %d\n",
	               pnum, jcol, w, nrow);
            #endif
            }
        };
    };
};
